cc sage from GB30A    at FLINDERS.EDU.AU      on 06/06/88 at 13:53:01

C   IMSL ROUTINE NAME   - CLEQT
C
C-----------------------------------------------------------------------
C
C   COMPUTER            - VAX/SINGLE
C
C   LATEST REVISION     - JUNE 1, 1981
C
C   PURPOSE             - MATRIX DECOMPOSITION, LINEAR EQUATION
C                           SOLUTION - SPACE ECONOMIZER SOLUTION -
C                           COMPLEX MATRICES
C
C   USAGE               - CALL CLEQT (A,N,IA,B,M,IB,IJOB,WA,IER)
C
C   ARGUMENTS    A      - INPUT COMPLEX N BY N MATRIX CONTAINING THE
C                           COMPLEX COEFFICIENTS OF THE EQUATION AX = B.
C                         ON OUTPUT, A CONTAINS THE L-U DECOMPOSITION OF
C                           A ROWWISE PERMUTATION OF THE INPUT MATRIX A.
C                N      - ORDER OF MATRIX A. (INPUT)
C                IA     - ROW DIMENSION OF MATRIX A EXACTLY AS
C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE
C                           CALLING PROGRAM. (INPUT)
C                B      - INPUT COMPLEX N BY M MATRIX CONTAINING THE
C                           M COMPLEX VALUED RIGHT HAND SIDES OF THE
C                           EQUATION AX = B.
C                         ON OUTPUT, THE SOLUTION MATRIX X REPLACES B.
C                           IF IJOB=1, B IS NOT USED.
C                M      - NUMBER OF RIGHT HAND SIDES (COLUMNS IN B).
C                           (INPUT)
C                IB     - ROW DIMENSION OF MATRIX B EXACTLY AS SPECIFIED
C                           IN THE DIMENSION STATEMENT IN THE CALLING
C                           PROGRAM. (INPUT)
C                IJOB   - INPUT OPTION PARAMETER.  IJOB=I IMPLIES WHEN
C                           I=0, FACTOR THE MATRIX AND SOLVE THE
C                             EQUATION AX=B.
C                           I=1, FACTOR THE MATRIX A.
C                           I=2, SOLVE THE EQUATION AX=B.  THIS
C                             OPTION IMPLIES THAT CLEQT HAS ALREADY
C                             BEEN CALLED USING IJOB=0 OR 1 SO THAT
C                             THE MATRIX HAS ALREADY BEEN FACTORED.  IN
C                             THIS CASE OUTPUT MATRIX A MUST HAVE BEEN
C                             SAVED FOR REUSE IN THE CALL TO CLEQT.
C                WA     - WORK AREA OF LENGTH N CONTAINING THE PIVOT
C                           INDICES.
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         TERMINAL ERROR
C                           IER=129 INDICATES THAT MATRIX A IS
C                             ALGORITHMICALLY SINGULAR.  (SEE THE
C                             CHAPTER L PRELUDE.)
C
C   PRECISION/HARDWARE  - SINGLE AND DOUBLE/H32
C                       - SINGLE/H36,H48,H60
C
C   REQD. IMSL ROUTINES - VERTST,VGETIO
C
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP
C
C   REMARKS  1.  WHEN IJOB=1, ARGUMENTS B, M AND IB ARE NOT USED BY
C                CLEQT.
C            2.  INPUT MATRIX A IS DESTROYED WHEN IJOB=0 OR 1. WHEN
C                IJOB=0 OR 2, B IS REPLACED WITH THE SOLUTION X.
C            3.  CLEQT CAN BE USED TO COMPUTE THE INVERSE OF A COMPLEX
C                MATRIX. THIS IS DONE BY CALLING CLEQT WITH M=N,
C                B=THE N BY N IDENTITY MATRIX AND IJOB=0. WHEN N IS
C                LARGE, IT MAY BE MORE PRACTICAL TO COMPUTE THE INVERSE
C                A COLUMN AT A TIME. TO DO THIS, FIRST CALL CLEQT WITH
C                IJOB=1 TO FACTOR A. MAKE SUCCEEDING CALLS WITH M=1, B
C                =A COLUMN OF THE IDENTITY MATRIX AND IJOB=2. B WILL BE
C                REPLACED BY THE CORRESPONDING COLUMN OF A INVERSE.
C            4.  THE DETERMINANT OF A CAN BE COMPUTED AFTER CLEQT HAS
C                BEEN CALLED AS FOLLOWS
C
C                  DET = (1.0,0.0)
C                  DO 5 I = 1,N
C                     IPVT = WA(I)
C                     IF (IPVT .NE. I) DET = -DET
C                     DET = DET*A(I,I)
C                5 CONTINUE
C
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.
C
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE CLEQT (A,N,IA,B,M,IB,IJOB,WA,IER)
C                                  SPECIFICATIONS FOR ARGUMENTS
      implicit real*8(a-h,o-z)
      INTEGER*4          N,IA,M,IB,IJOB,IER
      COMPLEX*16            A(IA,N),B(IB,M)
      REAL*8             WA(N)
C
C                               SPECIFICATIONS FOR LOCAL VARIABLES
C
      REAL*8             P,Q,ZERO,ONE,T(2),RN,BIG
      DOUBLE COMPLEX     SUM
      COMPLEX*16            TEMP,TEMP2
      INTEGER            I,J,JM1,IM1,K,IMAX,JP1,IW,N1
      EQUIVALENCE        (TEMP2,T(1))
      DATA               ZERO/0.0E00/,ONE/1.0E0/
C
C                                  INITIALIZATION
C                                  FIRST EXECUTABLE STATEMENT
C
      IER = 0
      IF (IJOB .EQ. 2) GO TO 75
      RN = REAL(N)
C
C                                  FIND EQUILIBRATION FACTORS
C
      DO 10 I = 1,N
         BIG = ZERO
         DO 5 J = 1,N
            TEMP = A(I,J)
            P = ABS(TEMP)
            IF (P .GT. BIG) BIG = P
    5    CONTINUE
         IF (BIG .EQ. ZERO) GO TO 105
         WA(I) = ONE/BIG
   10 CONTINUE
C
C                  L-U DECOMPOSITION
C
      DO 70 J = 1,N
         JM1 = J - 1
         IF (JM1 .LT. 1) GO TO 25
C
C                  COMPUTE U(I,J), I=1,...,J-1
C
         DO 20 I = 1,JM1
            SUM = A(I,J)
            IM1 = I - 1
            IF (IM1 .LT. 1) GO TO 20
            DO 15 K = 1,IM1
               SUM = SUM - A(I,K)*A(K,J)
   15       CONTINUE
            A(I,J) = SUM
   20    CONTINUE
   25    P = ZERO
C
C                                  COMPUTE U(J,J) AND L(I,J), I=J+1,...,
C
         DO 45 I = J,N
            SUM = A(I,J)
            IF (JM1 .LT. 1) GO TO 40
            DO 35 K = 1,JM1
               SUM = SUM - A(I,K)*A(K,J)
   35       CONTINUE
            A(I,J) = SUM
   40       Q = WA(I)*ABS(SUM)
            IF (P .GE. Q) GO TO 45
            P = Q
            IMAX = I
   45    CONTINUE
C
C                                  TEST FOR ALGORITHMIC SINGULARITY
C
         Q = RN + P
         IF (Q .EQ. RN) GO TO 105
         IF (J .EQ. IMAX) GO TO 60
C
C                                  INTERCHANGE ROWS J AND IMAX
C
         DO 50 K = 1,N
            TEMP = A(IMAX,K)
            A(IMAX,K) = A(J,K)
            A(J,K) = TEMP
   50    CONTINUE
         WA(IMAX) = WA(J)
   60    WA(J) = IMAX
         JP1 = J + 1
         IF (JP1 .GT. N) GO TO 70
C
C                                  DIVIDE BY PIVOT ELEMENT U(J,J)
C
         TEMP = A(J,J)
         DO 65 I = JP1,N
            A(I,J) = A(I,J)/TEMP
   65    CONTINUE
   70 CONTINUE
   75 IF (IJOB .EQ. 1) GO TO 9005
      DO 103 K = 1,M
C
C                                  BACKSUBSTITUTION
C                                  SOLVE UX = Y FOR X
C
         IW = 0
         DO 90 I = 1,N
            IMAX = WA(I)
            SUM = B(IMAX,K)
            B(IMAX,K) = B(I,K)
            IF (IW .EQ. 0) GO TO 85
            IM1 = I - 1
            DO 80 J = IW,IM1
               SUM = SUM - A(I,J)*B(J,K)
   80       CONTINUE
            GO TO 88
   85       CONTINUE
            TEMP2 = SUM
            IF (T(1) .NE. ZERO .OR. T(2) .NE. ZERO) IW = I
   88       CONTINUE
            TEMP2 = SUM
            B(I,K) = TEMP2
   90    CONTINUE
C
C                                  SOLVE LY = B FOR Y
C
         N1 = N + 1
         DO 100 IW = 1,N
            I = N1 - IW
            JP1 = I + 1
            SUM = B(I,K)
            IF (JP1 .GT. N) GO TO 98
            DO 95 J = JP1,N
               SUM = SUM - A(I,J)*B(J,K)
   95       CONTINUE
   98       B(I,K) = SUM/A(I,I)
  100    CONTINUE
  103 CONTINUE
      GO TO 9005
C
C                                  ALGORITHMIC SINGULARITY
C
  105 IER = 129
 9000 CONTINUE
C                                  PRINT ERROR
      CALL VERTST(IER,6HCLEQT )
 9005 RETURN
      END
C   IMSL ROUTINE NAME   - ITCLEQ
C
C-----------------------------------------------------------------------
C
C   COMPUTER            - VAX/SINGLE
C
C   LATEST REVISION     - NOVEMBER 1, 1984
C
C   PURPOSE             - LINEAR EQUATION SOLUTION - COMPLEX MATRIX -
C                           HIGH ACCURACY SOLUTION
C
C   USAGE               - CALL ITCLEQ (A,N,IA,B,M,IB,IJOB,WA,WK,IER)
C
C   ARGUMENTS    A      - INPUT COMPLEX MATRIX OF DIMENSION N BY N
C                           CONTAINING THE COMPLEX COEFFICIENTS OF THE
C                           EQUATION AX = B.
C                N      - ORDER OF MATRIX A. (INPUT)
C                IA     - ROW DIMENSION OF MATRIX A EXACTLY AS
C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE
C                           CALLING PROGRAM. (INPUT)
C                B      - INPUT COMPLEX MATRIX OF DIMENSION N BY M
C                           CONTAINING THE M COMPLEX VALUED RIGHT HAND
C                           SIDES OF THE EQUATION AX = B.
C                         ON OUTPUT, THE SOLUTION MATRIX X REPLACES B.
C                           IF IJOB = 1, B IS NOT USED.
C                M      - NUMBER OF RIGHT HAND SIDES (COLUMNS IN B).
C                           (INPUT)
C                IB     - ROW DIMENSION OF MATRIX B EXACTLY AS
C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE
C                           CALLING PROGRAM. (INPUT)
C                IJOB   - INPUT OPTION PARAMETER. IJOB = I IMPLIES,
C                           I = 0, FACTOR THE MATRIX AND SOLVE THE
C                             EQUATION AX = B.
C                           I = 1, FACTOR THE MATRIX A.
C                           I = 2, SOLVE THE EQUATION AX = B. THIS
C                             OPTION IMPLIES THAT ITCLEQ HAS ALREADY
C                             BEEN CALLED USING IJOB = 0 OR 1 SO THAT
C                             THE MATRIX HAS ALREADY BEEN FACTORED. IN
C                             THIS CASE WORK AREAS WA AND WK MUST HAVE
C                             BEEN SAVED FOR REUSE IN THE CALL TO
C                             ITCLEQ.
C                WA     - COMPLEX WORK AREA OF LENGTH N*(N+2). THE
C                           FIRST N*N LOCATIONS CONTAIN THE LU
C                           DECOMPOSITION OF A ROWWISE PERMUTATION OF
C                           A. THE REMAINDER OF WA IS USED AS WORK
C                           SPACE.
C                WK     - WORK AREA OF LENGTH N. ON OUTPUT WK CONTAINS
C                           THE PIVOT INDICES.
C                IER    - ERROR PARAMETER. (OUTPUT)
C                         TERMINAL ERROR
C                           IER=129 INDICATES THAT MATRIX A IS
C                             ALGORITHMICALLY SINGULAR. (SEE THE
C                             CHAPTER L PRELUDE).
C                           IER=130 INDICATES THAT ITERATIVE
C                             IMPROVEMENT FAILED TO CONVERGE. THE
C                             MATRIX IS TOO ILL CONDITIONED.
C
C   PRECISION/HARDWARE  - SINGLE AND DOUBLE/H32
C                       - SINGLE/H36,H48,H60
C
C   REQD. IMSL ROUTINES - SINGLE/CLEQT,VERTST,VGETIO
C                       - DOUBLE/CLEQT,VERTST,VGETIO,VXADD,VXMUL,
C                           VXSTO
C
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP
C
C   REMARKS  1.  WHEN IJOB=1, ARGUMENTS B, M AND IB ARE NOT USED BY
C                ITCLEQ.
C            2.  WHEN IJOB=0 OR 2, B IS REPLACED WITH THE SOLUTION X.
C            3.  ITCLEQ CAN BE USED TO COMPUTE THE INVERSE OF A COMPLEX
C                MATRIX. THIS IS DONE BY CALLING ITCLEQ WITH M=N,
C                B=THE N BY N IDENTITY MATRIX AND IJOB=0. WHEN N IS
C                LARGE, IT MAY BE MORE PRACTICAL TO COMPUTE THE INVERSE
C                A COLUMN AT A TIME. TO DO THIS, FIRST CALL ITCLEQ WITH
C                IJOB=1 TO FACTOR A. MAKE SUCCEEDING CALLS WITH M=1, B
C                =A COLUMN OF THE IDENTITY MATRIX AND IJOB=2. B WILL BE
C                REPLACED BY THE CORRESPONDING COLUMN OF A INVERSE.
C            4.  THE DETERMINANT OF A CAN BE COMPUTED AFTER ITCLEQ HAS
C                BEEN CALLED AS FOLLOWS
C
C                  DET = (1.0,0.0)
C                  DO 5 I = 1,N
C                     IPVT = WK(I)
C                     IF (IPVT .NE. I) DET = -DET
C                     INDX = I + (I-1)*N
C                     DET = DET*WA(INDX)
C                5 CONTINUE
C
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.
C
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE ITCLEQ  (A,N,IA,B,M,IB,IJOB,WA,WK,IER)
C
      implicit real*8(a-h,o-z)
      COMPLEX*16            A(IA,N),B(IB,M),WA(N,1),TEMPA,TEMPB
      REAL*8             WK(N),TA(2),TB(2)
      REAL*8             AR,AI,BR,BI,CR,CI,DXNORM,XNORM,ZERO
      REAL*8             ANORM,BNORM
      DOUBLE PRECISION   DACC
      EQUIVALENCE        (TA(1),TEMPA),(TB(1),TEMPB),
     *                   (TA(1),AR),(TA(2),AI),(TB(1),BR),(TB(2),BI)
      DATA               ZERO/0.0E0/
      DATA               ITMAX/40/
C
C                        FIRST EXECUTABLE STATEMENT
C
      IER = 0
      N1 = N + 1
      N2 = N + 2
      IF (IJOB .EQ. 2) GO TO 15
C
C     SAVE MATRIX A, COMPUTE THE NORM OF A AND B.
C
      ANORM = ZERO
      BNORM = ZERO
      DO 10 I = 1,N
      TEMPB = B(I,1)
      BNORM = AMAX1(BNORM,ABS(BR),ABS(BI))
         DO 5 J = 1,N
            WA(I,J) = A(I,J)
            TEMPA = A(I,J)
            ANORM = AMAX1(ANORM,ABS(AR),ABS(AI))
    5    CONTINUE
   10 CONTINUE
C
C          FACTOR MATRIX A
C
      CALL CLEQT (WA,N,N,B,M,IB,1,WK,IER)
      IF (IER .NE. 0) GO TO 9000
      IF (IJOB .EQ. 1) GO TO 9005
C
C                                  SAVE THE RIGHT HAND SIDES
C
   15 CONTINUE
      DO 65 J = 1,M
         DO 20 I = 1,N
            WA(I,N1) = B(I,J)
   20    CONTINUE
C
C                OBTAIN AN INITIAL SOLUTION
C
         CALL CLEQT(WA,N,N,WA(1,N1),1,N,2,WK,IER)
C
C                COMPUTE THE NORM OF THE SOLUTION
C
         XNORM = ZERO
         DO 25 I = 1,N
            TEMPA = WA(I,N1)
            XNORM = AMAX1(XNORM,ABS(AR),ABS(AI))
   25    CONTINUE
         IF (XNORM .EQ. ZERO) GO TO 65
C
C       START OF THE ITERATIVE REFINEMENT LOOP
C       FIRST COMPUTE RESIDUALS
C
         DO 50 ITER = 1,ITMAX
            DO 40 I = 1,N
               TEMPB = B(I,J)
               DACC = BR
               DO 30 JJ = 1,N
                  TEMPA = A(I,JJ)
                  TEMPB = WA(JJ,N1)
                  DACC = DACC - AR*BR
                  DACC = DACC + AI*BI
   30          CONTINUE
               CR = DACC
               TEMPB = B(I,J)
               DACC = BI
               DO 35 JJ = 1,N
                  TEMPA = A(I,JJ)
                  TEMPB = WA(JJ,N1)
                  DACC = DACC - AR*BI
                  DACC = DACC - AI*BR
   35          CONTINUE
               CI = DACC
               WA(I,N2) = CMPLX(CR,CI)
   40       CONTINUE
            CALL CLEQT(WA,N,N,WA(1,N2),1,N,2,WK,IER)
            DXNORM = ZERO
C
C     UPDATE THE SOLUTION AND TEST FOR CONVERGENCE BY COMPUTING
C     THE NORM OF THE DIFFERENCE BETWEEN THE TWO SOLUTIONS.
C
            DO 45 I = 1,N
               WA(I,N1) = WA(I,N1) + WA(I,N2)
               TEMPA = WA(I,N2)
               DXNORM = AMAX1(DXNORM,ABS(AR),ABS(AI))
   45       CONTINUE
            IF (XNORM+DXNORM .EQ. XNORM) GO TO 55
   50    CONTINUE
         IER = 130
C
C             STORE THE SOLUTION
C
   55  CONTINUE
         DO 60 JK = 1,N
            B(JK,J) = WA(JK,N1)
   60    CONTINUE
         IF (IER .NE. 0) GO TO 9000
   65 CONTINUE
      GO TO 9005
 9000 CONTINUE
      CALL VERTST(IER,6HITCLEQ )
 9005 RETURN
      END
C
C
C   PROGRAM NAME : VERTST (IMSL ROUTINE MODIFICATION)
C
*   VERSION      : VAX/SINGLE                                        *
*                                                                      *
*   DATE         : 06 JUNE 1986                                        *
*                                                                      *
*   USAGE        : CALL VERTST (IER,NAME)                              *
*                                                                      *
*   ARGUMENTS    :  INPUT(S) : IER : ERROR PARAMETER:                  *
*                                    IER =  I+J WHERE I = 128 IMPLIES  *
*                                    TERMINAL  ERROR MESSAGE; I =  64  *
*                                    IMPLIES  WARNING  WITH FIX MESS-  *
*                                    AGE;   I =  32  IMPLIES  WARNING  *
*                                    MESSAGE. J = ERROR CODE RELEVANT  *
*                                    TO CALLING ROUTINE.               *
*                             NAME : A CHARACTER STRING OF LENGTH SIX  *
*                                    PROVIDING THE  NAME OF THE CALL-  *
*                                    ING ROUTINE.                      *
*                                                                      *
*   CALL(S) TO   : VGETIO,VSPKD                                        *
*                                                                      *
*                                                                      *
*                            NOTES ON USAGE                            *
*                            --------------                            *
*                                                                      *
*   PURPOSE - PRINT A MESSAGE REFLECTING AN ERROR CONDITION.           *
*   NOTATION  -  INFORMATION  ON SPECIAL NOTATION AND CONVENTIONS IS   *
*   AVAILABLE IN THE MANUAL INTRODUCTION.                              *
*                                                                      *
*   REMARKS      THE ERROR  MESSAGE PRODUCED BY  VERTST IS WRITTEN TO  *
*                THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT NUMBER CAN  *
*                BE DETERMINED BY CALLING VGETIO AS FOLLOWS:           *
*                                                                      *
*                            CALL VGETIO (1,NIN,NOUT).                 *
*                                                                      *
*                THE  OUTPUT  UNIT  NUMBER  CAN BE CHANGED BY CALLING  *
*                VGETIO AS FOLLOWS:                                    *
*                                                                      *
*                            NIN = 0                                   *
*                            NOUT = NEW OUTPUT UNIT NUMBER             *
*                            CALL VGETIO(3,NIN,NOUT)                   *
*                                                                      *
*                SEE THE VGETIO DOCUMENT FOR MORE DETAILS.             *
*                                                                      *
*   COPYRIGHT           - 1982 BY IMSL, INC. ALL RIGHTS RESERVED.      *
*                                                                      *
************************************************************************
*
      SUBROUTINE VERTST (IER,NAME)
*
*  SPECIFICATIONS FOR ARGUMENTS
*
      INTEGER            IER
      INTEGER            NAME(1)
*
*  SPECIFICATIONS FOR LOCAL VARIABLES
*
      INTEGER            I,IEQ,IEQDF,IOUNIT,LEVEL,LEVOLD,NAMEQ(6),
     :                   NAMSET(6),NAMUPK(6),NIN,NMTB
*
      DATA               NAMSET/1HU,1HE,1HR,1HS,1HE,1HT/
      DATA               NAMEQ/6*1H /
      DATA               LEVEL/4/,IEQDF/0/,IEQ/1H=/
*
*  UNPACK NAME INTO NAMUPK
*
*  FIRST EXECUTABLE STATEMENT
*
      CALL VSPKD (NAME,6,NAMUPK,NMTB)
*
*  GET OUTPUT UNIT NUMBER
*
      CALL VGETIO (1,NIN,IOUNIT)
*
*  CHECK IER
*
      IF (IER .GT. 999) GO TO 25
      IF (IER .LT. -32) GO TO 55
      IF (IER .LE. 128) GO TO 5
      IF (LEVEL .LT. 1) GO TO 30
*
*  PRINT TERMINAL MESSAGE
*
      IF (IEQDF .EQ. 1) WRITE (IOUNIT,35) IER,NAMEQ,IEQ,NAMUPK
      IF (IEQDF .EQ. 0) WRITE (IOUNIT,35) IER,NAMUPK
      GO TO 30
    5 IF (IER .LE. 64) GO TO 10
      IF (LEVEL .LT. 2) GO TO 30
*
*  PRINT WARNING WITH FIX MESSAGE
*
      IF (IEQDF .EQ. 1) WRITE (IOUNIT,40) IER,NAMEQ,IEQ,NAMUPK
      IF (IEQDF .EQ. 0) WRITE (IOUNIT,40) IER,NAMUPK
      GO TO 30
   10 IF (IER .LE. 32) GO TO 15
*
*  PRINT WARNING MESSAGE
*
      IF (LEVEL .LT. 3) GO TO 30
      IF (IEQDF .EQ. 1) WRITE (IOUNIT,45) IER,NAMEQ,IEQ,NAMUPK
      IF (IEQDF .EQ. 0) WRITE (IOUNIT,45) IER,NAMUPK
      GO TO 30
   15 CONTINUE
*
*  CHECK FOR UERSET CALL
*
      DO 20 I=1,6
         IF (NAMUPK(I) .NE. NAMSET(I)) GO TO 25
   20 CONTINUE
      LEVOLD = LEVEL
      LEVEL = IER
      IER = LEVOLD
      IF (LEVEL .LT. 0) LEVEL = 4
      IF (LEVEL .GT. 4) LEVEL = 4
      GO TO 30
   25 CONTINUE
      IF (LEVEL .LT. 4) GO TO 30
*
*  PRINT NON-DEFINED MESSAGE
*
      IF (IEQDF .EQ. 1) WRITE (IOUNIT,50) IER,NAMEQ,IEQ,NAMUPK
      IF (IEQDF .EQ. 0) WRITE (IOUNIT,50) IER,NAMUPK
   30 IEQDF = 0
      RETURN
   35 FORMAT(19H *** TERMINAL ERROR,10X,7H(IER = ,I3,
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)
   40 FORMAT(27H *** WARNING WITH FIX ERROR,2X,7H(IER = ,I3,
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)
   45 FORMAT(18H *** WARNING ERROR,11X,7H(IER = ,I3,
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)
   50 FORMAT(20H *** UNDEFINED ERROR,9X,7H(IER = ,I5,
     1       20H) FROM IMSL ROUTINE ,6A1,A1,6A1)
*
*  SAVE P FOR P = R CASE
*  P IS THE PAGE NAMUPK
*  R IS THE ROUTINE NAMUPK
*
   55 IEQDF = 1
      DO 60 I=1,6
   60 NAMEQ(I) = NAMUPK(I)
   65 RETURN
*
      END
C
C
*   PROGRAM NAME : VGETIO (IMSL ROUTINE MODIFICATION)                  *
*                                                                      *
*   VERSION      : VAX/SINGLE                                        *
*                                                                      *
*   DATE         : 09 JUNE 1986                                        *
*                                                                      *
*   USAGE        : CALL VGETIO (IOPT,NIN,NOUT)                         *
*                                                                      *
*   ARGUMENTS    :  INPUT(S) :IOPT : INPUT OPTION                      *
*                                                                      *
************************************************************************
*
      SUBROUTINE VGETIO (IOPT,NIN,NOUT)
*
*  SPECIFICATIONS FOR ARGUMENTS
*
      INTEGER            IOPT,NIN,NOUT
*
*  SPECIFICATIONS FOR LOCAL VARIABLES
*
      INTEGER            NIND,NOUTD
*
      DATA               NIND/5/,
     :                   NOUTD/1/
*
*  FIRST EXECUTABLE STATEMENT
*
      IF (IOPT .EQ. 3) GO TO 10
      IF (IOPT .EQ. 2) GO TO 5
      IF (IOPT .NE. 1) GO TO 9005
      NIN = NIND
      NOUT = NOUTD
      GO TO 9005
    5 NIND = NIN
      GO TO 9005
   10 NOUTD = NOUT
 9005 RETURN
*
      END
C
C
*   PROGRAM NAME : VSPKD (IMSL ROUTINE MODIFICATION)                   *
*                                                                      *
*   VERSION      : VAX/SINGLE                                        *
*                                                                      *
*   DATE         : 04 JUNE 1986                                        *
*                                                                      *
*   USAGE        : CALL VSPKD (PACKED,NCHARS,UNPAKD,NCHMTB)            *
*                                                                      *
*   ARGUMENT(S)  :  INPUT(S) :PACKED: CHARACTER STRING TO BE UNPACKED. *
*                             NCHARS: LENGTH OF PACKED. SEE REMARKS.   *
*                                                                      *
*                  OUTPUT(S) :UNPAKD: INTEGER  ARRAY  TO  RECEIVE THE  *
*                                     UNPACKED REPRESENTATION  OF THE  *
*                                     STRING.                          *
*                             NCHMTB: NCHARS MINUS TRAILING BLANKS.    *
*                                                                      *
*                           NOTES ON USAGE                             *
*                           --------------                             *
*                                                                      *
*   NUCLEUS  CALLED  BY  IMSL  ROUTINES  THAT  HAVE  CHARACTER STRING  *
*   ARGUMENTS                                                          *
*                                                                      *
*   VSPKD UNPACKS A CHARACTER STRING INTO AN INTEGER ARRAY IN (A1)     *
*   FORMAT. UP TO 129  CHARACTERS MAY BE USED.  ANY IN EXCESS OF THAT  *
*   ARE IGNORED.                                                       *
*                                                                      *
*   COPYRIGHT           - 1984 BY IMSL, INC.  ALL RIGHTS RESERVED.     *
*                                                                      *
************************************************************************
*
      SUBROUTINE VSPKD  (PACKED,NCHARS,UNPAKD,NCHMTB)
*
*  SPECIFICATIONS FOR ARGUMENTS
*
      INTEGER            NC,NCHARS,NCHMTB,UNPAKD(1),IBLANK,PACKED(1)
*
      DATA               IBLANK /1H /
*
*  INITIALIZE NCHMTB
*
      NCHMTB = 0
*
*  RETURN IF NCHARS IS LE ZERO
*
      IF(NCHARS .LE. 0) RETURN
*
*  SET NC=NUMBER OF CHARS TO BE DECODED
*
      NC = MIN0 (129,NCHARS)
cx    DECODE (NC,150,PACKED) (UNPAKD(I),I=1,NC)
  150 FORMAT (129A1)
*
*  CHECK UNPAKD ARRAY AND SET NCHMTB BASED ON TRAILING BLANKS FOUND
*
      DO 200 N = 1,NC
         NN = NC - N + 1
         IF(UNPAKD(NN) .NE. IBLANK) GO TO 210
  200 CONTINUE
      NN = 0
  210 NCHMTB = NN
*
      RETURN
      END
